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A description of models and boundary conditions required for coupling radiation and 
ablation physics to a hypersonic flow simulation is provided. Chemical equilibrium routines 
for varying elemental mass fraction are required in the flow solver to integrate with the 
equilibrium chemistry assumption employed in the ablation models. The capability also 
enables an equilibrium catalytic wall boundary condition in the non-ablating case. The 
paper focuses on numerical implementation issues using FIRE II, Mars return, and Apollo 
4 applications to provide context for discussion. Variable relaxation factors applied to the 
Jacobian elements of partial equilibrium relations required for convergence are defined. 
Challenges of strong radiation coupling in a shock capturing algorithm are addressed. 
Results are presented to show how the current suite of models responds to a wide variety 
of conditions involving coupled radiation and ablation. 


I. Nomenclature 


Bold face, lowercase variable names refer to vectors. Bold face, uppercase variable names refer to matri- 
ces. Bracketed entry indicates units or quantity used to non-dimensionalize. 


Roman symbols 

A-j ~ Ej 

Ci,abl 

Dj 

f 

frad 

F 

F 

x 1,1 
AF° 

Gj k 

H n 

Kj{T) 


L 


ref 


m 


curve fit constants for Kj(T) to define equilibrium relation for non-base species j 
elemental mass fraction of element i in ablation products 
effective diffusion coefficient for species j, [L re f\ G>] 
species flux vector 

fraction of divergence of radiative flux introduced into flow solver (Eq 19) 
transformation matrix from species continuity to elemental continuity equations 
element of F defined in Eq. 1 
change in free energy in Eqs. 7 and 8, [J/kmole] 

coefficient of In pk in partial equilibrium relation for non-base species j (Eq. 10) 

enthalpy of species n, [J /kmole] 

equilibrium constant for reaction j, Eq. 11 

reference length used to non-dimensionalize distance [m] 

blowing rate of ablation products per unit area, [pooRoo] 
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Mi 

Mj 

N exchange 
N rad 

Pn 

P 

q 

q 

R 

R 

S n 

t 

T 

w 

x,y 
u, v 
Vi 

V, 

V , 30 
z 


molecular weight of element i, [kg/kmole] 
molecular weight of species j, [kg/kmole] 
frequency for updating ghost cells 

frequency for updating divergence of radiative flux in interior 
partial pressure of species n, [atm] 
pressure, [pooV^] 

vector of dependent variables, species densities [poo] 
heating rate, [W/m 2 ] 

universal gas constant, 8314.3 [ J/K -kmole] 
modified gas constant defined below Eq. 9 
entropy of species n, [J /kmole-K] 
time [f^re//koo] 
temperature [K] 

vector of chemical source terms, [pooVoo / L re f] 
distance, Cartesian coordinates [L re f\ 
velocities in x and y directions, respectively, [V^o] 
diffusion velocity of element i, [V/o] 
diffusion velocity of species j, [V/o] 
reference velocity in free stream [m/s] 

10000. /T used in Eq. 11 


Greek symbols 

a i,j 

Pm,j 

Pn,j 

e 

£ 

Pi 

Pj 

Poo 

a 

Xj 


number of atoms of element i in species j 

stoichiometric coefficient of reactant m in non-base species relation j 
stoichiometric coefficient of product n in non-base species relation j 
surface emissivity 

damping coefficient on surface temperature updates 

density of element i, [Poo] 

density of species j, [poo] 

reference density in free stream, [kg/m 3 ] 

Stefan-Boltzmann constant, 5.66961 10 -8 W/(m 2 -K 4 ) 
mole fraction of species j 


Subscripts 

0 

1 

i 

j 

k 

TO 

n 

oo 


ghost cell behind boundary 

internal cell adjacent to boundary 

element index 

species index 

species index 

reactant species index 

product species index 

reference condition in free stream 


II. Introduction 

The design of many planetary entry vehicles requires definition of the aerothermal environment including 
shock-layer radiation and ablation. In air, shock-layer radiation becomes significant at velocities above 10 
km/s for vehicles of order 1 meter diameter or greater. Radiative heating is proportional to the shock-layer 
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thickness (to some positive fractional power less than or equal to 1 depending on optical thickness of the 
gas), which increases as the vehicle diameter increases. The radiation (photon transport) cools the shock 
layer but also directs energy toward the vehicle surface to increase the heating rate. In environments where 
radiative heating is significant, an ablating thermal protection system is usually required to protect the 
vehicle (unless some active cooling system is employed). The ablation process (sublimation and reactions 
involving solid char and pyrolysis of other resins composing the virgin material) divert a significant fraction 
of the convective and radiative energy flux away from the heat shield bond line. Examples of missions that 
require(d) some level of coupled radiation and ablation analyses include Apollo, Orion, Stardust, Pioneer 
Venus, and Galileo. 

The numerical simulation of coupled radiation and ablation has progressed over the past 30 years from 
viscous shock layer methods 1 to upwind based finite-volume methods. 2-4 Implementation details for the 
present work on coupled radiation and ablation differ in some important ways. In recent papers describing 
coupled radiation and ablation, Matsuyama 2 simplifies the system by modeling only two effective species, 
gas from the free stream and ablation gas which are each characterized by constant elemental mass fraction. 
Fujita 3,4 includes a comprehensive thermochemical nonequilibrium flow model in conjunction with a species 
mass continuity boundary condition that includes reaction probabilities at the wall and a specified pyrolysis 
gas species mass fraction blowing rate. The present method assumes chemical equilibrium at the wall with 
the net elemental flux from blowing and diffusion balanced by the specified elemental mass fraction of 
the ablation products. An option is retained in which gas chemistry in the shock layer can be calculated at 
equilibrium by employing partial equilibrium relations to supplement elemental continuity or by full chemical 
nonequilibrium employing species continuity. In either case, the original species densities are retained as 
fundamental dependent variables. 

This paper is a companion to the recent papers by Johnston et. al. 5,6 Whereas these companion papers 
focus on specific applications, the present paper focuses on the modeling and the numerical details of model 
implementation. The baseline flow solver is a newly refactored version of the LAURA code' that utilizes 
Fortran 95 modules. The physics modules described herein (as well as the refactored LAURA) are a subset 
of the FUN3D 8 suite of modules so that they are available to that unstructured flow solver. 

III. Free Energy Minimization - (FEM) 


A. Overview 

The Gibbs Free Energy minimization (FEM) option is implemented with a goal of making simple adjustments 
to the baseline thermochemical nonequilibrium flow solver to enable a thermochemical equilibrium solver 
for generic mixtures of gases with spatially varying elemental mass fractions. The capability is required to 
enable interfaces with ablation routines based on equilibrium gas chemistry. The capability also enables 
simulations that are so deep into the equilibrium regime that the magnitude of species chemical source 
terms overwhelm the magnitude of convection and diffusion. (Roundoff error in the evaluation of the source 
terms is significant relative to the magnitude of convective and diffusive terms.) Consequently, the sum of 
the species continuity equations does not telescope to a meaningful mixture continuity equation because of 
roundoff error. The modification is implemented as follows. 

1. A parser is added to extract the constituent elements from the species names. 

2. Basis species are identified as the simplest homogeneous particle (smallest molecular weight, usually 
the atom) corresponding to each constituent element. 

3. Partial equilibrium relationships (derived from free energy minimization) for each remaining species 
are based on the following criteria. 

(a) Neutral molecules are equilibrated to constituent base species 

(b) Ionized heavy particles are equilibrated with electron and corresponding neutral base. 

4. The species continuity equations corresponding to base species in the baseline thermochemical nonequi- 
librium solver are replaced with elemental conservation equations. The elemental conservation equa- 
tions are formed from an appropriate linear combination of species conservation equations already 
present in the baseline algorithm. The electron continuity equation is replaced by an algebraic relation 
among species to force net charge neutrality. 
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5. The species continuity equations corresponding to non-base species in the baseline thermochemical 
nonequilibrium solver are replaced with algebraic, partial equilibrium relations from step 3. 

6. The baseline solver (point-implicit or line-implicit relaxation) is applied to the modified system of 
equations. Even though elemental conservation equations are present the dependent variable set retains 
the use of the original set of species, never using elemental densities directly. 

The algorithm described here has limitations. The algorithm fails if any element does not appear as a 
homogeneous component of one of the species. The algorithm does not handle solid carbon in the mix of 
species although an equilibrium state with solid carbon at the surface can be utilized to compute sublimation 
temperature. The partial equilibrium relations potentially can be very stiff (large magnitude of transformed 
equilibrium constant K :) defined below) because no consideration is given to local conditions that may govern 
a better selection of equilibrating species. In fact, a floor on equilibrium constant is applied which kicks in 
at temperatures less than 1000 K. The impact of this floor is that the mass fractions of trace species may 
be in error at low temperatures but the consequence of the error is in the noise. The algorithm could be 
modified to address any of these limitations but in tests conducted to date the simple approach described 
here is sufficient. 


B. Elemental Continuity Equations 

The transformation from species continuity to elemental continuity involves the identification of the fraction 
of element i in each species j. The transformation matrix F involves N e rows, one for each element, and 
N s columns, one for each species. The fraction is computed by multiplying the elemental weight, M, by the 
number of atoms of that element in the species, ctij and dividing by the species weight M j . 


_ dtjj Mi 

M,- 

Consider the one-dimensional species conservation equations set. 

<9q <9f 

1 = w 

dt dx 

where f includes convection and diffusion. If one multiplies both sides of Eq. 
conservation equations 


(1) 

(2) 

2 by F one obtains the elemental 


F 


<9q 

dt 


or 

dx 


<9q di 

dt + dx 


Fw 


0 


( 3 ) 

( 4 ) 


Note that Fw = 0. Chemical reactions cannot create or destroy an element in the control volume. Still, it is 
important to apply the transformation matrix before adding the source term because the critical components 
of continuity through convection and diffusion may be of the order of roundoff error in the assembly of the 
source term w. The assembly of the Jacobian of the elemental conservation laws (Eq. 12) employs the same 
matrix transform and retains the use of species densities (rather than elemental densities) as the fundamental 
dependent variables. Consequently, the remaining non-base species must be replaced by a partial equilibrium 
relation in order to close the equation set for species densities. 


C. Partial Equilibrium Relations 

The derivation of partial equilibrium relations through the minimization of free energy is available in several 
texts. 9 Consider the partial equilibrium relation for non-base species j which may be expressed as 


reactants products 

Y Pm,jM m < = > Y Pn,jN n 

m n 


( 5 ) 
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where one of the products N n corresponds to non-base species j, the reactants correspond to constituent base 
species, and one additional species may be required to enforce a chemical balance. The partial equilibrium 
relation involving non-base species j may then be derived as follows. 




RT 

n 

a S n | , RT ^ ^ a ( 

/ jPn,j{ jjrj, + l n M ) 


AF j V- 0 , RT ^ 0 , RT 


M n 

m 

Hm - TS„ 
RT 


m^) 

M ' 


Kj{T) 


-RT ln 

I^T n— products P n 

(6) 

_ Wm— reactants P 171 

^ ^ fin,j In Pn y ^ fim,j In p m (7) 

n 

m 


y ^ fin,j ln p n + y ^ fim,j ln Pm (8) 

n 

species 

m 


y Gyfclnpfc 

(9) 


k 


where R = pooR/ 101325. to include conversion from pressure in atmospheres to the metric system and to 
non-dimensionalize density by poo- The non-base species continuity equations are thus replaced by a linear 
equation in ln p^ with coefficients defined 


! 0 if species k is not part of reaction j 

— fin • j if species k = n is a product in reaction j (10) 

Pm,j if species k = m is a reactant in reaction j 

The temperature dependent equilibrium constant Kj(T) is curve fit using 

Kj (T) = A~ + Bj + Cj \nZ+ DyZ + EjZ 2 (11) 

Z 

where Z = 10000. /T and the constants A 3 - Ej are defined by fitting to the left-hand-side of Eq. 8 at 
T = 1000., 2000., 4000., 8000., 16000. The enthalpy H n and entropy S n for species n is obtained from curve 
fits. 10 

The coupled solution of Eq. 9 with the other conservation laws can be a challenge, particularly when it 
is early in the simulation and the solution is far from a converged state. It has been found that ignoring the 
Jacobian of Eq. 9 with respect to temperature and by multiplying the Jacobian of Eq. 9 with respect to 

species densities by a safety factor of 15 has enabled solution of the FEM coupled equation set for all tests 

conducted for this paper. The factor 15 was determined through numerical experiment. 


IV. Ablation Boundary Conditions 

The ablation boundary condition derived here assumes there will be coupling with an ablation module 
such as FIAT 11 in conjunction with an equilibrium chemistry module. It is assumed that ablation mass flow 
rate m = rh pyro i ys i s + m c har , the elemental composition of ablation gases Cj >a b;, and the temperature of the 
ablating surface are supplied by the ablation module or can be user specified. The ablation module itself 
requires the aerothermodynamic environments (pressure, radiative and convective heating, and possibly 
shear) over time to compute material response. Consequently, a boundary condition must be applied to 
the flow solver that closes the equation sets for either nonequilibrium external flow using species continuity 
equations or equilibrium external flow using elemental continuity equations. As noted in the previous section 
the species densities are retained as the fundamental dependent variable data set regardless of the assumption 
of equilibrium state of the external flow. 

It is assumed that gas chemistry at and within the ablating heatshield is in an equilibrium state. A 
steady, one-dimensional, elemental continuity equation (Eq. 12 ) describes the flow of an equilibrium gas 
across the boundary that accounts for convection and diffusion of element i. 

dpjV _ dpiVj 

dy dy ( j 
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Because we wish to work directly with species densities the transformation matrix introduced in the previous 
section is applied in Eq. 13 below. 


dpjv 

J dy 


= F, 


*,3 


9pjVj 

dy 


(13) 


where summation over species j is implied. 

Integrate Eq. 13 over y from some virtual point sufficiently deep in the heatshield where it is assumed 
that the net elemental mass flux is given by Ci^bWT- This assumption ignores the possibility that ablating 
mass flux sources are distributed through the depth of the char. Short of simulating flow through the char 
it is not clear that this modeling assumption can be improved. The resulting relation (Eq. 14) provides a 
balance for each element in which the convected elemental mass flux out of the boundary minus the diffused 
elemental mass flux into the boundary must equal the ablated elemental mass flux specified for the material. 


FijPjV -Ci, a birh = FjjpjVj (14) 

■ 

F i,j PjV Ci : ablril — F ijpDj ( 1 ^) 

Note that an effective binary diffusion coefficient is applied to model diffusive flux. A mass corrected diffusive 
flux 12 is also applied to guarantee that the sum of diffusive flux over all species (or over all elements) goes 
to zero. Simulations show some elemental separation near the wall for the no blowing case. Elemental 
separation in the no-blowing case is less when modeling multi-component diffusion using the Stefan Maxwell 
equations. 12 A verification check shows that no elemental separation is computed in the no-blowing case 
when a constant Schmidt number is applied. 

The discretization of Eq. 15 is shown in Eq. 16. It is necessary to solve for all dependent variables in 
the ghost cell 0. 


'(Pjv) l + (pjv) o 




(-i.ablFl — Fj,'j 


(pF>j) i + (pDj) o 


Xj, i Xj,o 

Ay 


(16) 


The solution in the ghost cell requires simultaneous relaxation of the partial equilibrium relations given by 
Eq. 9. The Jacobian formed using In pj as the dependent variable provides a more robust convergence of this 
boundary condition. A relaxation factor is applied to each row of the Jacobian according to the magnitude of 
the current residual for that row equal to n/ (dln[pj} 2 + n) where n is the local iteration count on relaxation 
steps for this boundary condition. Note that in the limit of zero blowing this boundary condition specifies 
an equilibrium catalytic wall - that is, the species are in chemical equilibrium as determined by the wall 
temperature and pressure and elemental mass fraction gradient equal to zero. 

\ [(Hi + (pu)o] = rh (17) 


Pi + pw\ =p 0 + Poi’o (18) 

The equation set is closed with a discretization of net mass flux (Eq. 17) and normal momentum 
balance (Eq. 18). Eq. 17 is used to eliminate Vo from all equations. Then Eq. 18 replaces one of the 
elemental continuity equations so that the ghost cell pressure po is consistent with the specified temperature 
To = 2 T w — Ti and the computed densities pjj). 

At present, ablation rates and ablation gas elemental mass fractions are explicitly set by the user. The 
boundary condition is sufficiently robust to enable a coupled simulation assuming equilibrium gas chemistry 
at Voa =15 km/s and a dimensionless blowing rate equal to 0.2. 

On the first step of initialization, all cells (including ghost cells) are set to inflow conditions. On the 
second step of initialization ghost cell velocities are reset to zero and wall temperatures are set to a user 
specified value. Simulations involving ablation and/or radiation are started from a converged (or nearly 
converged) non-radiating and non-ablating solution. The line-implicit formulation for relaxation of interior 
cells requires a Jacobian of the boundary condition. It is assumed for this purpose that species mass fractions 
and temperature in the ghost cells are fixed. (The point-implicit formulation for relaxation of interior points, 
by definition, requires no Jacobian information of neighbor cells.) 
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During interior cell updates the surface convective heating, shear, and diffusion flux are calculated from 
current values in the ghost cells. Ghost cell values in turn are updated at a user defined frequency, N exc hange- 
(Ghost cells are also used to hold information from neighbor blocks and the variable name N exc hange derives 
from the frequency in which neighbor block information is exchanged.) Ghost cell values of species density 
are fully defined by interior point values - there is no time dependent update or damping using any fraction 
of a previous value. Damping is applied to the update of a radiative equilibrium surface temperature in the 
case of no blowing. The radiative equilibrium wall temperature boundary condition is then defined with 
T w new = £[(7/(ecr)] 1 / 4 + (1 — e)T w old . In this case, the ghost cell temperature is computed as T 0 = 2 T w — T\. 

Nexchange is set equal to 2 in all cases presented herein. The temperature damping factor e is varied 
between 0.001 and 0.01. In more recent applications (some involving a surface energy balance not covered 
here) the combination N exc h an ge = 20 and e = 0.1 have provided more robust convergence. 

V. Radiation Coupling 

Radiative energy transport is computed using the tangent slab approximation with the HARA code. 13,14 
The flow solver provides species densities, temperatures, and distances along a ray orthogonal to the surface 
as inputs to HARA. HARA returns divergence of the radiative flux and radiative flux to the wall for use in 
LAURA. Calls are made to HARA on the order of hundreds of relaxation steps, N rad , holding the radiative 
quantities constant over this period (loosely coupled). Options are provided to interpolate radiative quantities 
between selected rays to reduce overall time in the radiation solver. The HARA radiation model treats atomic 
lines and continuum of N, O, C, and H. Molecular bands are treated for N 2 ,N^, 02,NO, C 2 , C3, CN, CO, C 2 H. 
Non-Boltzmann electronic state models are applied for atoms and molecules. Typical computational time 
for a single line-of-sight is 20 to 50 seconds - short enough to enable fully coupled simulations with codes 
like LAURA. 

A parameter f rad may be specified to introduce the divergence of the radiative flux more slowly into the 
flow solver. 


v • q ^ 1 = fradN • q r ad (HARA) + (1 - f rad )X7 ■ q" od (19) 

If the radiation introduces strong cooling into the shock layer the shock standoff distance decreases 
significantly. As the captured shock moves closer to the body, computational cells previously in the hot shock 
layer are overtaken by the pre-shock domain. However, these cells retain the cooling term V • q ra d,i,j,k from 
the previous call to the radiation solver. The LAURA code enforces a floor on temperature for every point in 
the flowfield that allows the simulation to survive such transients under most circumstances. However, cases 
have been observed where introduction of the full update to radiation has caused stability problems. Best 
practice so far has been to set f ra d to 0.4 on the first call from an uncoupled solution and then to increase 
its value to 1 on subsequent returns from the radiation solver. 

VI. Results and Discussion 

Three test cases are presented to demonstrate that the algorithms for coupled multiphysics (radiation plus 
ablation plus flow solver) indeed work and that they provide credible simulations (consistent with existing 
flight data and benchmark simulations). All cases were run on an Intel based Mac laptop, each simulation 
requiring on the order of two to four hours for 10000 to 20000 relaxation steps. Courant numbers varied 
from 5 to 50 (local time stepping). (More recent work suggests that constant time stepping with time steps 
equal to time for free stream to travel 0.001 m is more efficient.) Further details of each case are provided 
in the original sources. 

A. Fire II, Coupled Radiation Only 

The Fire II data set from 1967 remains the touchstone for validation of entry heating simulations of radiative 
heating. 15-17 The vehicle had three overlayed beryllium heatshields. A fresh surface was exposed twice during 
entry after the outermost heatshield was jettisoned as its integrity became compromised due to thermal load. 
Consequently, radiation data was obtained without the complications of ablation products. Metallic surfaces 
are generally considered to be strongly catalytic but the potential formation of an oxide coating introduces 
uncertainty in the proper formulation of this boundary condition. Therefore results are presented using 
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bounding approximations for a super-catalytic condition (species mass fractions on the surface equal to 
free stream values) and for a non-catalytic condition (zero species fraction gradient at the wall). Surface 
temperatures were specified according to flight data. 



Front 


Side 


Figure 1. Sensors on Fire II vehicle, circle represents radiometer locations, x represents calorimeter locations. 


Instrument locations are presented in Fig. 1. The radiometers record radiation intensity contribution 
between 0.4 and 6.2 eV (window limits). The calorimeters record total convective and absorbed radiative 
heating. The absorbed radiative energy is calculated as a function of wavelength. 18 Simulations at two 
trajectory points are intended to illustrate capabilities at a strong thermochemical nonequilibrium condition 
(1636 seconds) and a near-equilibrium condition (1643 seconds). 

The computed heating rates for the subject points are compared to experimental data in Fig. 2. The 
total heating rates measured by the calorimeters are presented in Figs. 2 (a) and (b). The solid black 
line in these figures indicate the uncoupled total heating (radiative heating computed from converged, non- 
radiating case) assuming a super-catalytic boundary condition. The solid red line indicates the total heating 
with full coupling of radiative energy transfer. The coupling reduces the total heating by 5 to 10 % for these 
trajectory points. The dashed lines indicate the coupled, non-catalytic heating rate. The catalytic effect 
is most prominent (reduction by 1/3) at the strong nonequilibrium point at 1636 seconds. The catalytic 
effect is still evident at 1643 seconds with a net reduction of 5 %. The experimental data, indicated by 
blue circles, is fully bounded by the extremes in modeling of surface catalysis. The uncoupled radiation 
intensity is presented with a dashed line and the coupled radiation intensity is presented with a solid red 
line in Figs. 2 (c) and (d). The uncertainty in the experimental data in these figures (blue circles and error 
bar) bounds the computed coupled results. The simulations indicate that radiation coupling is required to 
match experimental data and that the suite of physical models used to compute these cases are consistent 
with experimental data. 

B. Mars Return, Coupled Radiation and Ablation, Equilibrium Shock Layer 

The Mars return example has been extensively studied for over 30 years, beginning with the then state-of- 
the-art viscous shock layer (VSL) code methodology. 1- 19,20 The cases are useful because shock layer profiles 
and parametric studies are available for comparison in the literature. The VSL simulations used a fitted 
shock boundary and equilibrium chemistry. The LAURA simulations are the first to work this case with a 
shock capturing algorithm. This case also demonstrates the Free Energy Minimization routine with elemental 
continuity equations replacing the species continuity equations. Free stream conditions for flow over a 3.05 
m radius sphere are Vx = 15.24 km/s, p x = 0.000255 kg/m 3 . Elemental mass fractions of injected gas are 
[C, H, O, N] = (0.92, 0.022, 0.049, 0.009). The wall temperature is set to 3600 K and the dimensionless 


8 of 12 


American Institute of Aeronautics and Astronautics Paper 2009-1399 



(a) Total heating distribution at 1636 seconds. 



(c) Radiation intensity distribution at 1636 seconds. 



(b) Total heating distribution at 1643 seconds. 



(d) Radiation intensity distribution at 1643 seconds. 


Figure 2. Total and radiative heating distributions around Fire II vehicle at trajectory points 1636 and 1643. 


blowing rate to is set to 0. (no ablation) and 0.2. 

Stagnation streamline profiles of temperature and wall directed radiative flux computed by LAURA / 
HARA (solid lines) and VSL 1 (dashed lines) are compared in Fig. 3. Three cases are compared: (1) non- 
ablating, non-radiating; (2) non-ablating, coupled radiating; and (3) coupled ablating, coupled radiating. 
The non-radiating profile shows a constant temperature across the inviscid portion of the shock layer, as 
expected for an equilibrium gas. The shock layer cooling with radiation and associated decrease in shock 
standoff distance is evident in the middle set of curves of Fig. 3 (a). Note that the dashed lines from Moss 
in this figure end abruptly at the fitted shock used in the VSL algorithm. The ablating case shows a large 
increase in the boundary layer thickness (ss 3 cm) and corresponding displacement of the bow shock. The 
wall directed radiative flux in Fig. 3 (b) indicates absorption of photon flux in the ablation products near 
the wall. In all cases the LAURA / HARA simulations are in good agreement with the earlier VSL results 
- demonstrating consistency with a previous formulation of this coupled physics simulation. More details of 
this case, including comparisons with other methods, are available in the literature. 5 

C. Apollo 4, Coupled Radiation and Ablation 

The Apollo 4 test case, recently reviewed by Park, 21 provides flight data on radiative heating including 
effects of ablation products. Free stream conditions are V x = 10.25 km/s and = 0.000341 kg/m 3 on a 
3. meter radius sphere. Elemental mass fractions of injected gas are [C, H, O, N] = (0.63, 0.06, 0.30, 0.01). 
These fractions were derived from data discussed by Park 21 in which quasi-steady ablation is assumed and 
the elemental fractions of char and pyrolysis gas are combined according to the fraction of their individual 
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(b) Wall directed radiative flux profile 


Figure 3. Stagnation streamline profiles for Mars return benchmark case. 


blowing rates. The wall temperature is set to 2500 K and the dimensionless blowing rate to is 0.0086. 21 
Thermochemical nonequilibrium model is applied. The case was initialized from a non-ablating solution in 
which the equilibrium catalytic wall boundary condition was applied. 




Figure 4. Surface heating distributions for Apollo 4 entry with and without radiation and ablation coupling. 


Fig. 5 (a) plots the cumulative intensity arriving at the radiometer starting from 0.4 eV and concluding 
at 6.2 eV (the radiometer window is opaque outside this range) . The path length indicated by the solid 
line includes the shock-layer thickness immediately above the radiometer. The path length indicated by the 
dashed line includes the shock-layer thickness plus 8 cm of cavity depth (Fig. 5 (b)) which is assumed to 
be filled with gas at conditions on the surface. (There was no flow field simulation that included a cavity.) 
From 0.4 to approximately 2.0 eV there is essentially no difference in the intensity reaching the radiometer 
over either path length. The intensity reaching the radiometer above 2eV becomes path length dependent, 
indicating that some of the radiation above 2eV is being absorbed by the gas in the cavity. The circular 
symbol with uncertainty bar indicates the measured cumulative intensity. It is convenient to put this symbol 
at the end of the window range at 6.2 eV but the data in fact represents the integrated intensity from 0.4 to 
6.2eV. This figure shows that for the assumed blowing conditions the predicted cumulative intensity at the 
surface exceeds the measure value by approximately 20 %. However, if the effect of absorption over a longer 
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(b) Extra path length through cavity 
to radiometer 


Figure 5. Comparisons with Apollo 4 radiometer data. 


path length through the cavity is included then excellent agreement with measurement is obtained. 

We caution that this agreement does not constitute validation. The specified blowing rates and elemental 
mass fraction of the ablation products are estimates 21 - not computed on any first principles within this 
simulation. Still, the results confirm that for a reasonable estimate of boundary conditions there is consistency 
with measured data - providing more supporting evidence that the coupled radiation and ablation algorithms 
are correctly implemented. Further details of this simulation with more consideration of parametric variation 
is available in the literature. 5 


VII. Concluding Remarks 

A description of models and boundary conditions required for coupling radiation and ablation physics to 
a hypersonic flow simulation finite- volume solver is provided. The focus herein is to present the algorithmic 
details (relaxation factors, update frequency, approximations to the Jacobians) that are required to maintain 
a stable, fully coupled simulation. Three application examples are then provided (FIRE II, Mars return, and 
Apollo 4 applications) to illustrate that: (1) - the new modules are correctly implemented; and, (2) - the 
results are consistent with either experimental data or prior, benchmark calculations. 

Chemical equilibrium routines for varying elemental mass fraction are required in the flow solver to 
integrate with the equilibrium chemistry assumption employed in the ablation models. This capability is 
provided through a strong coupling of partial equilibrium relations derived through Free Energy Minimization 
(FEM). The FEM capability also enables an equilibrium catalytic wall boundary condition in the non- 
ablating case. 

The HARA code, using a tangent slab approximation, forms the basis of the radiation solver. Its relatively 
quick solution time per line of sight (20-50 seconds) makes the fully coupled solution process tenable. The 
ablation is currently modeled using a user specified mass flow rate and user specified elemental mass fractions 
of the blown gas. Work is ongoing to compute a char ablation rate based on equilibrium of gaseous carbon 
at the wall with solid carbon in the char and to couple this system with the FIAT 11 code to compute the 
pyrolysis gas rates. 

The FIRE II test cases include the entire forebody surface distributions and compare them to experimental 
heating data. The predicted heating rates for super-catalytic and non-catalytic boundary conditions bound 
the experimental data. The predictions for radiative intensity are within the uncertainty bounds of the 
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measured data. The Mars return test cases demonstrate the application of FEM over the entire shock layer 
including ablation products. Ablation coupling has the strongest impact on convective heating and milder 
impact on radiative heating in this case - indicative of the effect of absorption by ablation gases. Radiation 
coupling has the strongest impact on radiative heating and milder impact on convective heating. Computed 
results are consistent with benchmarks solutions obtained with a Viscous Shock Layer (VSL) approach. The 
Apollo 4 test case shows good agreement with a cumulative radiation intensity measurement through a shock 
layer including absorption through a boundary layer with ablation gases. These comparisons are collected 
here in one place to demonstrate code capability and consistency with available data for a challenging array 
of problems. 

Future work will engage an equilibrium state between atomic carbon and char along with a surface energy 
balance to compute a char blowing rate. Then coupling with FIAT will be introduced to handle the physics 
of pyrolysis gas without filtering through a blowing correction. A capability to automatically run a sequence 
of trajectory points from a user supplied file is currently running. It is planned to capture the heating history 
for input to FIAT. Finally, a quasi-one-dimensional grid adaptation capability will be modified to enable 
shape change as a function of the distributed blowing rates and the trajectory history. 
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